R5R Tutorial

library(ggplot2)

Introduction

In this tutorial you will be introduced to r5r: an interface to the R5 routing engine developed by Conveyal. R5 (r5r) allows you to do “rapid, realistic” routing on multimodal networks. It’s very flexible, enabling intermodal analyses and has the key strength of retaining the full detail of the public transport timetable. If you’re more familiar with python, you can consider using r5py instead - although it may not have all of the same functionality.

As an example, we’re going to be analyzing accessibility by public transport to beer gardens 🚇🍻

The tutorial consists of the following main steps

1. Installing r5r
2. GTFS preparation
3. Loading spatial data into R
4. Creating an Isochrone

In addition to r5r, you will be exposed to several other useful R packages: dplyr (for working with data tables), sf (for working with geodata), and tidytransit (for working with GTFS data). These are all very powerful packages that can make your life a lot easier if you learn how to use them! You’re encouraged to look into the documentation of these packages for more information.

Installing r5r

All the packages listed below (including r5r) can be installed by clicking “packages > install” on the sidebar and searching for the package name or by executing the following command: install.packages("package_name_here")

Since R5 (the routing engine itself) is written in java, we need to install the Java Development Kit to make r5r work. Specifically we need version 21. You can install it here: https://www.oracle.com/de/java/technologies/downloads/#java21

With that done, we can load all the packages we’ll be using today.

library(tidyverse)

Warning: package ‘tidyverse’ was built under R version 4.1.3Warning: package ‘tibble’ was built under R version 4.1.3Warning: package ‘tidyr’ was built under R version 4.1.3Warning: package ‘forcats’ was built under R version 4.1.3-- Attaching core tidyverse packages ------------------------------------------------------------ tidyverse 2.0.0 --
v dplyr     1.1.4     v readr     2.1.4
v forcats   1.0.0     v stringr   1.5.1
v ggplot2   3.5.1     v tibble    3.2.1
v lubridate 1.9.4     v tidyr     1.3.0
v purrr     1.0.4     -- Conflicts ------------------------------------------------------------------------------ tidyverse_conflicts() --
x dplyr::filter() masks stats::filter()
x dplyr::lag()    masks stats::lag()
i Use the ]8;;http://conflicted.r-lib.org/conflicted package]8;; to force all conflicts to become errors

library(sf)

Linking to GEOS 3.10.2, GDAL 3.4.1, PROJ 7.2.1; sf_use_s2() is TRUE

library(tidytransit)

options(java.parameters = '-Xmx4G') # we need to allocate RAM for r5r to work. You can increase this depending on your system resources.
library(r5r)

Please make sure you have already allocated some memory to Java by running:
  options(java.parameters = '-Xmx2G').
You should replace '2G' by the amount of memory you'll require. Currently, Java memory is set to -Xmx4G

GTFS Preparation

1. Download data from: https://gtfs.de/ - pick “Deutschland gesamt”
2. Move the zip file into the “raw_data” folder
3. Using the full feed for all of Germany in r5r will require a lot of resources. Since we’ll be analyzing a single city - lets filter the feed so that it is smaller. We’ll be using tidytransit for this.
4. Let’s first load in geodata representing the area we’ll be trimming to (in this case - Munich). We’ll use st_read from the sf package for this. It can handle all typical geospatial file types (.shp,.gpkg,etc.)

area <- st_read("./raw_data/munich_admin.gpkg")

Reading layer `boundary_administrative_munich' from data source 
  `C:\Users\barte\Documents\R\ilutm-r5r-tutorial\raw_data\munich_admin.gpkg' using driver `GPKG'
Simple feature collection with 25 features and 2 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 11.36078 ymin: 48.06162 xmax: 11.72291 ymax: 48.24812
Geodetic CRS:  WGS 84

5. It’s good practice to buffer our study area to avoid any edge effects. Let’s buffer our area by 5 km. NOTE: in this analysis we’re just going to be looking at accessibility to beer gardens in Munich. We aren’t considering any that are in the municipalities around the city. If your destinations extend beyond the border of your city, you may need to buffer your study area significantly! (>>5 km)
6. We’ve seen how to do this in QGIS, let’s use sf this time

area_buffered <- st_transform(area,crs = 25832) # switch crs to projected (ETRS89 / UTM zone 32N)

area_buffered <- st_union(area_buffered) # join geometries 

area_buffered <- st_buffer(area_buffered, 5000) # buffer by 5km

area_buffered <- st_transform(area_buffered,crs = 4326) #switch back to WGS84 for filtering the GTFS feed


ggplot()+
  geom_sf(data = area_buffered)+
  geom_sf(data = area)

7.Let’s load in the GTFS feed we downloaded and trim it. We can see that the number of stops in the feed has been significantly reduced.

8. Let’s save the filtered feed. R5 will use this to build a network, you can also use it for other public transport analyses (refer to previous tutorials). Tip: tidytransit has a lot of great functionality for this!

write_gtfs(gtfs_filtered,"./r5r_model/gtfs.zip")

r5r Set Up

1. Before running any more code we need to supply the correct data. So far we have prepared the GTFS data, but we also need a street network. R5 uses an extract of OSM data to build the street network. To do this, it needs to be provided with a .pbf file. This can be found on Geofabrik or https://slice.openstreetmap.us/. OSM by the slice is recommended because it has a cool name and because you can select the exact area you’re interested in (reducing network size, improving performance).
2. Just as when we were trimming the GTFS feed, extract an area larger than the city itself.

3. Save the .pbf file to the /r5r_model/ folder. The folder should look like this:

4. We can now tell r5 to build the network. Tip: if you are sure you aren’t modifying the GTFS or .pbf data, you can turn overwrite = F to avoid rebuilding the model (faster).

r5r_core <- setup_r5(data_path = "./r5r_model", verbose = FALSE,overwrite = T)

No raster .tif files found. Using elevation = 'NONE'.

Finished building network.dat at ./r5r_model/network.dat

OD Data

We need some data to use for origins and destinations. Our destinations will be beer gardens. The data is provided for you. It was prepared using an Overpass API query (as in previous tutorials).

beergardens <- st_read("raw_data/biergarten.gpkg")

Reading layer `centroids' from data source `C:\Users\barte\Documents\R\ilutm-r5r-tutorial\raw_data\biergarten.gpkg' using driver `GPKG'
Simple feature collection with 100 features and 3 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 11.3911 ymin: 48.07886 xmax: 11.71046 ymax: 48.21359
Geodetic CRS:  WGS 84

ggplot()+
  geom_sf(data = area)+
  geom_sf(data = beergardens)

For our origins, we’ll use public transport stops. As we discussed in lecture, a single “stop” may be represented by multiple in a GTFS feed (e.g., a bus stop on both sides of the street). We can take some steps to simplify this. We will group all stops that share the same value for parent_station and then keep a single coordinate by taking the average. Be aware that solving this in a line of code is often too good to be true. GTFS data often contains mistakes, and may require more careful cleaning!

stops <- summarize(group_by(gtfs_filtered$stops, parent_station),
                   stop_lat = mean(stop_lat),
                   stop_lon = mean(stop_lon),
                   parent_station_name = first(stop_name))  
stops

As a final step lets convert the table to geodata and filter it down to our study area.

stops <- st_as_sf(stops,coords = c("stop_lon","stop_lat"),crs = 4326,remove = F)

stops <-st_filter(stops,area) #apply spatial filter

stops

Simple feature collection with 1125 features and 4 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 11.38896 ymin: 48.06913 xmax: 11.71527 ymax: 48.22
Geodetic CRS:  WGS 84

ggplot()+
  geom_sf(data = area)+
  geom_sf(data = stops, color = "blue")+
  geom_sf(data = beergardens, size = 2)

beergardens

Simple feature collection with 100 features and 4 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 11.3911 ymin: 48.07886 xmax: 11.71046 ymax: 48.21359
Geodetic CRS:  WGS 84
First 10 features:
       osm_id    amenity                    name                      geom         id
1   110892783 biergarten Lindengarten Biergarten POINT (11.67838 48.11507)  110892783
2   368977873 biergarten         Tivoli Pavillon POINT (11.60341 48.16053)  368977873
3   502190647 biergarten             Olympia-Alm POINT (11.55432 48.17114)  502190647
4   292910957 biergarten       Tambosi Hofgarten  POINT (11.5786 48.14293)  292910957
5  1371548072 biergarten            Zum Biereder POINT (11.54993 48.16336) 1371548072
6    72273510 biergarten                Hirschau POINT (11.60248 48.16209)   72273510
7   343096349 biergarten          Paulanergarten  POINT (11.5718 48.20859)  343096349
8   325104410 biergarten   Ristorante Dolce Vita POINT (11.52478 48.08895)  325104410
9  2907073039 biergarten         Zur Gartenlaube POINT (11.60551 48.10014) 2907073039
10 1329567304 biergarten           La Locandiera  POINT (11.6409 48.08744) 1329567304

–> Before we start using r5r, lets format our data. The origins and destinations should have an id column. stops already has this, but beergardens does not. Let’s add a new column called id, and set the values to equal osm_id.

beergardens <- mutate(beergardens, id = osm_id)
stops <- mutate(stops, id = parent_station)

beergardens

Simple feature collection with 100 features and 4 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 11.3911 ymin: 48.07886 xmax: 11.71046 ymax: 48.21359
Geodetic CRS:  WGS 84
First 10 features:
       osm_id    amenity                    name                      geom         id
1   110892783 biergarten Lindengarten Biergarten POINT (11.67838 48.11507)  110892783
2   368977873 biergarten         Tivoli Pavillon POINT (11.60341 48.16053)  368977873
3   502190647 biergarten             Olympia-Alm POINT (11.55432 48.17114)  502190647
4   292910957 biergarten       Tambosi Hofgarten  POINT (11.5786 48.14293)  292910957
5  1371548072 biergarten            Zum Biereder POINT (11.54993 48.16336) 1371548072
6    72273510 biergarten                Hirschau POINT (11.60248 48.16209)   72273510
7   343096349 biergarten          Paulanergarten  POINT (11.5718 48.20859)  343096349
8   325104410 biergarten   Ristorante Dolce Vita POINT (11.52478 48.08895)  325104410
9  2907073039 biergarten         Zur Gartenlaube POINT (11.60551 48.10014) 2907073039
10 1329567304 biergarten           La Locandiera  POINT (11.6409 48.08744) 1329567304

TTM<-
travel_time_matrix(
  r5r_core,
  origins = stops, 
  destinations = stops, 
  mode = c("TRANSIT"),  
  mode_egress = "WALK", 
  departure_datetime = as.POSIXct("26-05-2025 9:00:00", format = "%d-%m-%Y %H:%M:%S"), #its important that the date and time is within your GTFS feed
  time_window = 10L, # 
  percentiles = c(10,50),
  max_walk_time = 10,
  max_trip_duration = 30L,
  walk_speed = 4,
  max_rides = 3,
)

TTM

Travel Time Matrix

The main functionality of r5r is the ability to calculate a travel time matrix. Knowing the travel time between pairs of points is useful on its own. It can also be used if you want more control over your accessibility calculations. For example, if you calculate a travel time matrix between pairs of public transport stops, you can calculate the access (origin –> first public transport stop) and egress (last public transport stop –> destination) legs of the trip independently. In other words, combining the output of r5r with the QGIS-based accessibility calculations we did in earlier tutorials.

Let’s give it a shot:

tip: there are many other settings you can change, you’re encouraged to read the documentation of r5r to see what is possible! The documentation has many nice examples!

• stops: our origins are the public transport stops
• destinations: our destinations are the public transport stops
• mode: the primary mode is public transport (which includes walking for access and transfers)
• mode_egress: the egress mode is walking
• departure_datetime: the time trips start
• time_window: set to 60. A travel time matrix will be calculated for each minute, for 60 minutes after the departure time. This allows for us to account for the variability in travel times. Imagine we ran our analysis only at 8:00 and we just missed a departure of an S-Bahn at 7:59, the travel time that we’d calculate would be severly penalized by excessive waiting. This parameter can be reduced (default is 10) as it can be computationally expensive. It’s increasingly important in suburban and rural areas where service frequencies are low.
• percentiles: this corresponds to the time_window. a value of 50 will report the median travel time. We can specify multiple percentiles by using a vector c()
• max_walk_time: max time for access, egress, and transfers by walking
• walk_speed: (kph)
• max_ride: maximum # of transfers

TTM <- mutate(TTM,tt_diff = travel_time_p50 - travel_time_p10)

print(paste0(round(sum(is.na(TTM$tt_diff)) / length(TTM$tt_diff)*100,0),"%"))

[1] "30%"

The impact of the time_frame can be quite significant! The percent of OD pairs that had a 50th percentile travel time > 30 min and a 10th percentile travel time <30 min:

TTM$tt_diff%>%na.omit()%>%summary()

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  0.000   2.000   3.000   3.127   4.000  21.000 

TTM<-TTM%>%select(-tt_diff)

Those that could be reached, were this much slower:

TTM

Simple feature collection with 294664 features and 5 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 11.38896 ymin: 48.06913 xmax: 11.71527 ymax: 48.22
Geodetic CRS:  WGS 84
First 10 features:
   from_id  to_id travel_time_p10 travel_time_p50 parent_station_name                  geometry
1   100060 100060               0               0   Josef-Lang-Straße POINT (11.47202 48.15011)
2   100060  10157              29              NA   Josef-Lang-Straße POINT (11.46016 48.18474)
3   100060 105596              28              NA   Josef-Lang-Straße POINT (11.48913 48.18399)
4   100060 109574              23              27   Josef-Lang-Straße POINT (11.47472 48.12805)
5   100060  11520              25              27   Josef-Lang-Straße POINT (11.54544 48.14428)
6   100060 128167              22              26   Josef-Lang-Straße POINT (11.46642 48.19177)
7   100060 131549              21              25   Josef-Lang-Straße POINT (11.44335 48.16361)
8   100060 131638              21              25   Josef-Lang-Straße POINT (11.47278 48.17485)
9   100060 132365              19              23   Josef-Lang-Straße  POINT (11.4531 48.14563)
10  100060   1335              27              NA   Josef-Lang-Straße POINT (11.45716 48.20257)

Let’s get the travel time matrix ready for an export. We will restructure the data so that the geometry of the destination stops (to_id) is added to each row. We will also add the name of the origin station so that the data is easier to work with. This is a similar operation to a VLOOKUP in excel. We use the from_id and to_id fields as they represent the parent_station.

For more guidance on dplyr functions (like left_join()) check out this cheat sheet

This format is easier to interpret if we take a single stop. For example, lets find all stops that can be reached from Theresienstraße within a median travel time of 20 minutes.

st_write(TTM,"./output/TTM.gpkg",append = F)

Writing layer `TTM' to data source `./output/TTM.gpkg' using driver `GPKG'
Writing 294664 features with 5 fields and geometry type Point.

Let’s export this so that we can work with it in QGIS.

st_write(TTM,"./output/TTM.gpkg",append = F)

We can achieve the same result as above by doing the following:

Then, by using what you’ve learned in previous tutorials, you can create service areas: